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1.  Introduction 


Since  the  work  in  the  late  60’s  by  B.  B.  Mandelbrot  and  J.  W.  Wallis  [1],  a 
novel  class  of  mathematical  functions  known  as  “fractals”  has  enjoyed 
great  popularity  [2-5].  Fractals  are  found  extensively  in  nature,  and  fractal 
surfaces  are  produced  by  a  number  of  basic  physical  processes  like  erosion, 
aggregation,  and  turbulent  flow.  These  processes  modify  shape  through  lo¬ 
cal  iterative  actions  and  thus  produce  fractals  as  an  end  result  [4,6,7]. 
Fractal  models  are  currently  in  use  in  a  number  of  areas  of  image  analysis, 
such  as  shape  from  texture  [7-9],  texture  analysis  and  segmentation 
[10,11],  interpolation  of  terrain  data  [12],  automatic  target  recognition 
[13,14],  and  medical  image  analysis  [15,26]  to  name  a  few.  The  realism  of 
computer-generated  images  using  fractal  models  strongly  suggests  that  the 
model  captures  relevant  perceptual  information  of  the  image  that  appeals  to 
our  visual  processing  system. 

1.1  Fractal  Model  and  Dimension 

A  fractal  can  be  loosely  defined  as  a  set  whose  only  consistent  description 
of  its  metric  properties  requires  a  noninteger  dimension  greater  than  or 
equal  to  the  standard  intuitive  definition  of  a  set’s  dimension.  This  leads  to 
the  notion  of  the  “fractal  dimension,”  which  is  the  consistent  measure  of  the 
set’s  roughness  and  irregularities  [7],  and  has  been  shown  to  be  twice  as  ac¬ 
curate  as  any  previous  reported  measure  for  perceptual  roughness  [17]. 
Fractals  are  classified  into  two  types:  deterministic  and  random,  and 
whereas  deterministic  fractals  show  self-similarity  over  all  scales,  random 
fractals  are  self-affine  over  a  limited  range  of  scales  [5].  A  sampling  of  the 
images  one  can  generate  with  either  type  can  be  found  in  the  references 
[2-5].  Random  fractals,  which  are  of  most  interest  to  the  scientific  commu¬ 
nity,  will  be  the  primary  concern  of  this  report. 

A  useful  mathematical  model  for  random  fractals,  as  proposed  by 
Mandelbrot  and  Van  Ness,  is  fractional  Brownian  motion  (fBm)  [1],  which 
is  a  generalization  of  the  familiar  Brownian  motion.  The  model  is  an  exten¬ 
sion  of  classical  Brownian  motion  and  encompasses  a  larger  group 
of  phenomena. 

Let  V(x)  be  a  self-affine /Bm  with  a  Euclidean  dimension  (E)  of  one,  with 
the  increments  having  zero-mean  Gaussian  distribution.  The  variance  has 
the  following  form: 


(|  V(x  +  £fr)  -  T(x)  f)  oe  I  iix:  1^^  ,  (1) 

and  has  a  fractal  dimension  D/ given  by  Df=E+l-H,  where  the  brackets 
“( ”  and  “ )”  denote  the  ensemble  averaged  over  many  samples  of  V(x).  H 
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is  often  called  the  “persistence”  parameter,  since  it  is  proportional  to  the  de¬ 
gree  of  correlation  and  is  bounded  above  by  one  and  below  by  zero.  The  in¬ 
crements  are  stationary  and  isotropic,  depending  only  on  the  difference  in  x, 
and  for  H  =  0.5  we  have  classical  Brownian  motion. 

(|V(j!:  +dr)- V(x)f)<x:|(5l[xj  .  (2) 

y(x)  shows  a  statistical  scaling  behavior;  i.e.,  if  x  is  changed  by  the  factor  k, 
then  the  increments  of  V(x)  change  by  the  factor  k  raised  to  the  H. 

{dV{kxf)  {dV{xf')k^^  .  (3) 

This  nonuniform  scaling  is  different  in  the  two  coordinates  x  and  y(x)  and 
gives  the_^m  its  self-affme  property,  unlike  the  self-similar  property  of  de¬ 
terministic  fractals  [27].  This  will  prove  important  in  the  selection  of  a  vi¬ 
able  fractal  dimension  estimator  [5]. 

The  extension  to  higher  dimensions  leads  to  a  fractal  Brownian  surface 
(fBs),  which  holds  our  current  interest  since  it  is  an  excellent  model  of  the 
earth’s  terrain  [7,12].  Let  y(x,y)  be  a  self-affine  JBs  with  a  Euclidean  di¬ 
mension  (E)  of  two  with  the  increments  having  zero-mean  Gaussian  distri¬ 
bution.  The  variance  again  has  the  following  form: 

(|  V{x  +  dx,y  +  dy)-  y(;c,y)  |^)  »=  |  dx^  +  dy^  \  ,  (4) 

and  has  a  fractal  dimension Z)/ given  by  Df=  E+\-H. 

The  characterization  of  a  generalized  power  spectrum  for  this  random  pro¬ 
cess  is  an  extension  of  the  concept  of  spectral  density  and  the  Weiner- 
Khintchine  relation  to  nonstationary  noise  [18-20]  and  takes  the  form  of  a 
random-phase  Fourier  transform,  with  the  expected  value  of  the  power 
obeying  the  relationship 


5(/)cc/-^  .  (5) 

Due  to  the  nonstationarity  of  fractal  Brownian  processes,  a  power  spectrum 
cannot  be  formally  defined;  however,  the  increments  (derivative)  of  the 
random  process  are  stationary,  and  therefore  a  power  spectrum  of  the  incre¬ 
ments  can  be  well-defined.  Using  this  property,  Flandrin  [21],  Keshner 
[22],  and  Womell  [23]  substantiated  the  above  form  of  the  power  spectrum. 
The  decay  parameter  B  is  related  to  the  fractal  dimension  and  to  the  persis¬ 
tence  parameter  H  by 


B  =2H  +  Ea.n<iDf  =  E+  1  -H, 


where  E  is  the  Euclidean  dimension  and  D/the  fractal  dimension.  If  any 
straight-line  cut  is  taken  in  the  x  or  y  direction  of  a  fractal  Brownian  sur¬ 
face,  the  spectral  density  will  have  the  power  law  in  equation  (5);  this  im¬ 
plies  that  a  two-dimensional  (2-D)  spectral  density  should  behave  as 
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/ 


2H  +  2 


(6) 


(  2 

for  spatial  frequencies  k  and  I,  f  =  \k  +l  )  . 

A  great  deal  of  work  has  been  done  in  both  synthesis  and  analysis  of 
fractals  using  this  power  spectrum  relationship  [9,13,28,29]. 


1.2  Synthetic  Fractal  Surfaces 

To  quote  H.  E.  Schepers  et  al  [24]  “production  of  a  truly  fractal  signal  is  it¬ 
self  a  nontrivial  exercise.”  This  truth  became  very  apparent  in  this  study  as 
I  found  that  the  only  reliable  measure  of  a  fractal  was  careful  development 
of  the  spectral  components,  comparison  to  other  authors’  images  and  sur¬ 
faces,  and  a  portion  of  faith.  The  two  algorithms  used  for  generating  fractal 
surfaces  were  the  mid-point  displacement  (MPD)  and  power  spectrum  den¬ 
sity  (PSD)  methods.  Both  are  currently  popular  in  the  computer  graphics 
community  for  synthetic  scene  generation  [5,16,27].  The  MPD  method  is  a 
recursive  technique  first  applied  by  Weiner  in  the  1920’ s  in  modeling 
Brownian  motion.  This  is  a  “quick  and  dirty”  algorithm  for  modeling 
fractals — quick  because  of  the  reduced  complexity  and  execution  speed  and 
dirty  since  it  sacrifices  mathematical  accuracy  (and  whose  result  can  be 
seen  as  visible  artifacts  on  a  surface  when  properly  shaded).  One  should  see 
Pietgen  and  Saupe  [5]  for  the  mechanics  of  this  algorithm  as  well  as  some 
interesting  renderings.  The  PSD  method  is  elaborated  on  by  Voss 
and  others  [5,27].  Basically,  one  creates  a  discrete  Fourier  transform  (DFT) 
of  the  real  image  u(m,n)  for  a  surface  with  the  2-D  DFT  pair  given  as 


V-l  N-\  .  , 

v{k,l)=  E  S  u{m,n)WNW!!," 

^  m=0  «=0 

0<k,l<N-l 

(7) 

,  N-l  A^-1  .  . 

u{m,n)  =  -^  X  E  v{k,l)WN'""m 
k=0  1=0 

0<m,  n<N-  1  , 

(8) 

(-2jn\ 
l%  =  expl  N  ) 

7 

(9) 

with  the  phase  information  derived  from  a  random  uniform  distribution  on 
27C.  When  the  DFT  coefficients  are  restricted  to  have  the  form 
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1 


(10) 


a  single  realization  of  the  fractal  surface  will  be  generated.  Alternatively, 
the  form 


H+\  ’ 


(11) 


will  generate  the  ensemble  average.  Here  H  =3-  Df  and  ( )  denote  expected 
values.  For  a  real  image  of  a  fractal  surface  we  must  also  observe  the  com¬ 
plex  conjugate  synunetry, 

v{k,l]  =  v*{N-k,N-l)  0<k,l<N-l  ;  (12) 

therefore,  our  phase  requirement — although  derived  from  a  random  uni¬ 
form  distribution  over  2k — has  to  be  zero  at  the  following  points:  (0,0), 
(0AV2),  {NI2,NI2),  and  (AV2,0)  [5,30]. 

One  can  readily  observe  that  this  model  imposes  a  bandlimited  spectrum, 
which  introduces  problems  in  producing  real  fractals  and  especially  those 
with  higher  fractal  dimension  values.  A  second  problem  is  the  periodic  na¬ 
ture  of  the  surface  image  created  by  the  inverse  fast  Fourier  transform 
(FFT).  This  periodicity  “perceptually  accentuates  the  presence  of  the  first 
spatial-frequency  harmonic”  [29]  of  the  DFT.  The  N/2  x  N/2  regions  of  the 
surface  were  used  for  the  analysis,  as  suggested  by  Voss  [31],  hopefully  to 
correct  this  problem. 

The  MPD  surfaces  will  serve  as  a  useful  comparison  to  the  more  reliable 
PSD  surfaces  in  both  the  fractal  dimension  calculations  and  three-dimen¬ 
sional  (3-D)  renderings  of  the  surface  images  and  DFTs.  Also,  one  must 
question  what  effect  the  artifacts  of  the  MPD  method  have  in  determining 
the  fractal  dimension  of  the  surfaces  created.  I  was  very  interested  in  the  re¬ 
sults  of  the  estimator  with  respect  to  the  single-realization  PSD  surfaces 
since  these  surfaces  are  more  representative  of  natural  relief  than  the  en¬ 
sembles  and,  therefore,  validate  the  estimator’s  performance. 


1.3  Fractal  Dimension  Estimation 

Several  techniques  have  been  developed  to  estimate  the  fractal  dimension 
as  researchers  probe  for  optimality  with  respect  to  a  particular  fractal  set 
[22,23,25].  All  the  techniques  attempt  to  measure  a  scaling  parameter  asso¬ 
ciated  with  either  a  self-similar,  statistically  self-similar,  or  self-affine 
fractal  (see  Pietgen  and  Saupe  [5]  for  more).  The  simplest  technique,  and 
quite  general  in  nature,  is  the  Mass  Scaling  (MS)  technique,  which  makes 
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no  assumptions  about  the  underlying  process  that  created  the  fractal.  In  the 
application  to  fractal  surfaces,  the  MS  technique  measures  the  mass  about 
some  central  point,  i.e.,  pixel  {x,y,z),  by  counting  the  number  of  pixels  M(l) 
that  fall  within  cubes  of  varying  side  1.  The  power  law  requirement  for  mass 
distribution  is  given  by 

M[l)  =  kff  ,  (13) 

where  £)/is  the  fractal  dimension  and  is  found  by  a  least-squares  fit  on  the 
plot  of  log  {M)  versus  log  (/)  [5,29,31]  and  k  is  the  proportionality  constant. 

Despite  improvements  to  the  earlier  work  of  Keller  et  al  [1 1],  the  MS  tech¬ 
nique  has  some  shortcomings;  first,  the  self-affine  property  creates  prob¬ 
lems  in  determining  the  fractal  dimension  since  scaling  may  be  different  in 
each  dimension  of  a  multidimension  space  as  elaborated  on  by  Voss  and 
others  [5,31].  Although  this  problem  is  cleverly  addressed  by  Moghaddam 
[29],  it  is  still  not  completely  resolved.  Secondly,  the  MS  method  is  not  re¬ 
liable  in  the  important  region-spanning  fractal  dimension  2.0  to  2.2,  which 
is  the  predominant  range  for  fractal  dimensions  of  the  earth’s  natural  relief 
[5,9]  and  is  completely  unreliable  at  very  high  fractal  dimensions  associ¬ 
ated  with  tree  clutter  [13,29]. 

The  power  spectrum  method  (PE)  is  based  on  the  spectral  model  of  fractal 
Brownian  motion,  which  is  more  complex  and  time  consuming,  but  has 
given  excellent  results  in  texture  segmentation  [7,28]  and  yields  the  least 
biased  results  and  lowest  variance  in  estimates  of  the  fractal  dimension  re¬ 
gardless  of  the  underlying  model  used  to  create  the  fractal  set  [24]. 

The  PE  estimator  in  this  report  will  determine  the  spectral  decay  parameter 
B  by  logarithmic  least-squares  fit  for  log  (power)  versus  log  (radial  fre¬ 
quency)  by  first  taking  a  local  8x8  window,  zero  padding  to  16  x  16,  and 
computing  the  EFT.  The  DFT  coefficients  are  used  to  calculate  the  total 
power  of  a  single  or  multiple  occurrence  of  a  specific  radial  frequency,  and 
then  a  logarithmic  plot  and  least-squares  fit  are  applied.  The  slope  (B)  is 
then  determined  and  used  to  calculate  the  fractal  dimension,  which  is 
mapped  back  to  the  central  four  pixels  of  the  8  x  8  window.  Finally,  the 
window  is  slid  by  two  pixels  and  the  process  repeated. 

The  type  of  windowing  function  applied  to  a  local  sample  of  the  test  surface 
is  an  area  of  debate  among  researchers.  The  Gabor  window  applied  to  the 
PE  method  has  proven  to  be  very  accurate  on  test  images  [9];  however,  this 
windowing  involves  the  extra  step  of  convolution  that  magnifies  execution 
time,  albeit  with  the  gain  of  improved  spatial-frequency  resolution  over  the 
domain  of  fractal  surfaces — i.e.,  surfaces  with  fractal  dimensions  between 
2.0  and  3.0.  In  this  study,  I  am  not  concerned  with  performance  in  the 
higher  fractal  dimension  regions  and,  in  reality,  optimality  is  not  measured 
solely  by  the  accuracy  of  the  algorithms  over  the  entire  range,  but  by  the 
practicality  of  implementation,  speed,  and  confidence  in  estimation  accu- 
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racy  over  limited  ranges.  Although  a  rectangular  window  has  shortcomings 
due  to  higher  sidelobe  energy,  the  implementation  is  simple  and  favors 
higher  frequency  resolution.  The  size  of  the  window  was  chosen  as  a  com¬ 
promise  between  spatial  resolution  improvement  over  larger  windows  and 
the  decreased  likelihood  of  including  nonstationary  regions  when  tested 
against  real  surfaces.  Also,  8x8  windows  seem  to  be  the  breakpoint  in  the 
variance  of  the  estimate,  with  windows  smaller  than  8x8  being  unreliable 
because  of  the  large  variance  of  the  estimated  fractal  dimension. 

The  zero  padding  effectively  gives  more  samples  of  the  2-D  DFT  envelope, 
which  is  useful  in  the  least-squares  fit  and  is  amenable  to  the  FFT  calcula¬ 
tion.  For  the  decay  parameter  B  calculation,  Pentland  [7]  performed  a  least- 
squares  fit  in  the  k  and  I  directions  of  the  discrete  spatial-frequency  plane 
only,  then  calculated  the  Bx  and  By  parameters,  followed  by  averaging 
these  B  values  to  determine  the  fractal  dimension.  He  relied  on  the  fact  that 
most  natural  fractals  exhibit  only  slight  anisotropy  [7,16].  Moghaddam  [29] 
improved  on  the  technique  by  using  a  weighted  least-squares  fit  based  on 
the  number  of  occurrences  of  a  particular  radial  frequency  in  the  discrete 
spatial-frequency  plane.  The  estimator  used  in  this  study  only  modified 
terms  due  to  a  single  occurrence  of  a  radial  frequency  (spatial  frequencies 
where  k  =  /),  giving  these  terms  equal  weight,  and  although  counter¬ 
intuitive,  gave  equivalent  results  as  a  simple  deweighting  algorithm  that 
removed  single  occurrenees  of  a  radial  frequency  in  the  DFT  before  the 
least-squares  fit.  The  mapping  to  the  central  four  pixels  is  based  on  the  ob¬ 
servations  of  Yokoya  et  al  [12].  Their  results  showed  that  although  the  in¬ 
tensity  of  the  fractal  surface  image  does  not  vary  smoothly,  the  fractal  di¬ 
mension  varies  smoothly  for  a  real  terrain  surface.  Therefore,  a  mapping  to 
the  central  four  pixels  should  not  decrease  the  reliablility  of  the  estimate 
significantly,  but  will  decrease  the  number  of  FFTs  performed  and  least- 
squares  calculations. 


2.  Procedure 


A  total  of  three  sets  of  fractal  surfaces  were  synthesized  using  (1)  the  MPD 
method,  (2)  the  PSD  method  with  the  DFT  magnitudes  derived  according  to 
the  ensemble  average  expression  (eq  (11)),  and,  finally,  (3)  the  PSD  method 
with  the  DFT  magnitudes  derived  according  to  the  single-realization  ex¬ 
pression  (eq  (10)).  The  fractal  dimensions  chosen  for  the  global  constructs 
were  2.0,  2.2,  2.4,  2.6,  2.8,  and  3.0.  The  same  seed  for  the  random  number 
generator  was  used  for  creating  all  synthetic  surfaces.  Global  constructs 
were  created  which  were  512x512,  then  segmented  into  four  blocks  of  256 
X  256  each.  This  led  to  a  total  of  72  surfaces  for  estimation  purposes. 

Estimation  of  the  fractal  dimension  was  performed  by  first  taking  a  local 
8x8  window,  zero  padding  to  16  x  16,  and  computing  the  FFT.  The  spec¬ 
tral  decay  parameter  B  was  determined  as  the  slope  of  the  logarithmic  least- 
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3.  Results 


Table  1.  Mean, 
variance,  and  mamse 
of  250  X  250  fractal 
dimension  arrays 
resulting  from  PE 
estimation  on  MPD 
constructs. 


squares  fit  for  log  (power)  versus  log  (radial  frequency)  and  then  used  to 
calculate  the  fractal  dimension,  which  is  mapped  back  to  the  central  four 
pixels  of  the  8  X  8  window.  Finally,  the  window  is  slid  by  two  pixels  and 
the  process  repeated.  This  generates  a  250  x  250  array  of  fractal  dimension 
values  for  each  256  x  256  fractal  surface.  First-  and  second-order  statistics 
were  calculated  for  each  fractal  dimension  array. 


A  total  of  72  surfaces  were  synthesized  and  subjected  to  tbe  PE  estimation 
algorithm.  Tables  1  to  3  list  the  three  classes  of  fractal  surfaces  created  and 
their  estimated  fractal  dimensions,  variance  of  the  estimate,  and  the  contri¬ 
bution  to  the  modified  average  mean  square  error  (mamse)  for  each  256  x 
256  block.  The  mamse  was  calculated  from  a  fractal  dimension  array  ac¬ 
cording  to  the  following  [30]: 

mamse  =47  S  \Df-Dfnf  ,  (14) 

/!=  1 

with  D/=  global  fractal  dimension  for  5 12  x  5 12  block  and  ^  =  the  mean 
fractal  dimension  for  256  x  256  subblock. 


File 

Fractal  dimension 
estimate 

Variance 

Modified 
average  mean 
square  error 

test201md.img 

2.3161 

0.0649 

0.0999 

test202md.img 

2.2135 

0.0775 

0.0456 

test203md.img 

2.0339 

0.0588 

0.0011 

test204md.img 

1.9885 

0.0493 

0.0001 

test221md.img 

2.4445 

0.0512 

0.0598 

test222md.img 

2.4440 

0.0681 

0.0595 

test223md.img 

2.2288 

0.0550 

0.0008 

test224md.img 

2.2161 

0.0499 

0.0003 

test241md.img 

2.6802 

0.0403 

0.0785 

test242md.img 

2.6783 

0.0542 

0.0774 

test243md.img 

2.5557 

0.0535 

0.0242 

test244mdamg 

2.5318 

0.0491 

0.0174 

test261md.img 

2.9331 

0.0337 

0.1109 

test262md.img 

2.9033 

0.0448 

0.0920 

test263md.img 

2.8769 

0.0481 

0.0767 

test264md.img 

2.8401 

0.0452 

0.0576 

test281md.img 

3.1536 

0.0304 

0.1250 

test282md.img 

3.1113 

0.0400 

0.0969 

test283md.img 

3.1251 

0.0413 

0.1057 

test284md.img 

3.1089 

0.0411 

0.0954 

test301md.img 

3.3195 

0.0287 

0.1021 

test302md.img 

3.2743 

0.0362 

0.0752 

test303md.img 

3.2945 

0.0359 

0.0867 

test304md.img 

3.2991 

0.0369 

0.0895 

Files  are  mid-point  displacement  realizations  with  fractal  dimensions 
2.0,  2.2,  2.4,  2.6,  2.8,  and  3.0. 
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Table  2.  Mean, 
variance,  and  mamse 
of  250  X  250  fractal 

File 

Fractal  dimension 
estimate 

Variance 

Modified 
average  mean 
square  error 

dimension  arrays 

n2psd201.img 

L9760 

0.0479 

0.0006 

resulting  from  PE 

n2psd202.img 

1.9529 

0.0272 

0.0022 

estimation  on  ensemble  n2psd203.img 

2.0356 

0.0275 

0.0013 

average  PSD 

n2psd204.img 

2.0107 

0.0498 

0.0001 

constructs. 

n2psd221.img 

2.0853 

0.0423 

0.1147 

n2psd222.img 

2.0377 

0.0240 

0.0263 

n2psd223.img 

2.1528 

0.0252 

0.0022 

n2psd224.img 

2.1479 

0.0430 

0.0027 

n2psd241.img 

2.3306 

0.0393 

0.0048 

n2psd242.img 

2.2766 

0.0239 

0.0152 

n2psd243.img 

2.3969 

0.0249 

0.0000 

n2psd244.img 

2,4086 

0.0385 

0.0001 

n2psd261.img 

2.6247 

0.0361 

0.0006 

n2psd262.img 

2.5894 

0,0242 

0.0001 

n2psd263.img 

2.6907 

0,0252 

0.0082 

n2psd264.img 

2.7121 

0.0358 

0.0126 

n2psd28Limg 

2.9198 

0.0343 

0,0144 

n2psd282.img 

2.9100 

0.0252 

0.0121 

n2psd283.img 

2.9794 

0.0260 

0.0322 

n2psd284.img 

3.0064 

0.0344 

0.0426 

/ 

n2psd301.img 

3.2008 

0.0348 

0.0403 

n2psd302.img 

3.2130 

0,0283 

0.0454 

n2psd303.img 

3.2408 

0.0284 

0.0580 

n2psd304,img 

3.2034 

0.0323 

0.0414 

Files  are  power  spectral  density  ensemble  average  realizations  with 

fractal  dimensions  2.0,  2.2,  2.4,  2.6,  2.8,  and  3.0. 

Table  3.  Mean, 

File 

Fractal  dimension 

Variance 

Modified 

variance,  and  mamse 

estimate 

average  mean 

of  250  X  250  fractal 

square  error 

dimension  array 

sr2psd201.img 

2.0740 

0.0562 

0.0055 

resulting  from  PE 

sr2psd202.img 

1.9762 

0.0322 

0.0006 

estimation  on  single 

sr2psd203.img 

2.1283 

0.0594 

0.0165 

realization  PSD 

sr2psd204.img 

2.0411 

0.0280 

0.0017 

constructs. 

sr2psd221.img 

2,1370 

0.0312 

0.0039 

sr2psd222.img 

2.2210 

0.0535 

0.0004 

sr2psd223,img 

2.2749 

0.0526 

0.0056 

sr2psd224.img 

2.2147 

0.0277 

0.0002 

sr2psd241.img 

2.3974 

0.0308 

0.0000 

sr2psd242.img 

2.4442 

0.0438 

0.0019 

sr2psd243,img 

2.4928 

0.0434 

0.0086 

sr2psd244.img 

2.4785 

0.0291 

0.0062 

sr2psd261.img 

2.6982 

0.0305 

0.0096 

sr2psd262,img 

2.7015 

0.0376 

0,0103 

sr2psd263.img 

2.7484 

0.0386 

0.0220 

sr2psd264.img 

2.7643 

0.0301 

0.0269 

sr2psd281.img 

2.9880 

0.0311 

0.0353 

sr2psd282.img 

2.9646 

0.0339 

0.0271 

sr2psd283.img 

3.0049 

0.0361 

0.0420 

sr2psd284.img 

3.0417 

0.0322 

0.0584 

sr2psd301.img 

3.2126 

0,0312 

0.0452 

sr2psd302.img 

3.2138 

0.0331 

0.0457 

sr2psd303.img 

3.2520 

0.0361 

0.0635 

sr2psd304.img 

3.2780 

0.0338 

0,0773 
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Files  are  power  spectral  density  single  realizations  with  fractal  dimensions 
2A2.2,  2.4,  2A  2.8,  and  3,0, 


The  trend  in  the  increase  of  the  mamse  in  both  the  single  realizations  and 
ensemble  averages  shows  the  difficulties  in  estimation  of  the  higher  dimen¬ 
sions  using  the  rectangular  window;  however,  the  results  at  lower  dimen¬ 
sions  are  very  encouraging.  The  low  variance  indicates  the  uniformity  of 
the  globally  generated  fractal  surfaces. 

Figure  1  plots  the  various  estimates  obteiined  against  fractal  dimension.  Fig¬ 
ures  2  to  5  are  representative  surface  plots  and  DFTs  of  the  single-realiza¬ 
tion  PSD  and  MPD  surfaces  produced  by  Gouraud  shading  with  an  over¬ 
head  light  source.  A  brief  remark  is  in  order  about  the  DFT  figures.  They 
were  created  by  taking  a  32  x  32  section  of  the  original  256  x  256  DFT 
(single  realization)  and  displaying  the  absolute  value  of  the  spatial  fre¬ 
quency  coefficients  from  {k,l)  =  (1,1)  to  ik,l)  =  (32,32);  the  axes  are  re¬ 
versed  since  the  sample  section  was  taken  from  the  DFT  region  (^,/)  = 
(224,224)  to  ikj)  =  (255,255),  which  is  analogous  by  symmetry. 

Closer  inspection  of  the  MPD  DFTs  suggests  that  extra  power  resides  along 
the  spatial  frequency  axes,  thus  affecting  the  PE  estimation  algorithm,  and, 
in  fact,  the  reliability  of  the  MPD  method  has  been  questioned  before  with 
regard  to  synthesizing  fractal  Brownian  motion  for  any  persistence  param¬ 
eter  H  deviating  from  1/2  [27].  The  slight  increase  in  the  power  residing 
along  the  spatial-frequency  axis  of  the  PSD  DFTs  can  be  accounted  for  by 
the  inherent  periodicity  of  the  construct;  therefore,  M4  x  NIA  samples  as 
used  in  Moghaddam’s  work  is  probably  preferred  and  may  have  led  to  a 
more  accurate  fractal  sample.  Figures  6  to  9  zire  representative  of  the  3-D 
surfaces  produced  by  the  different  synthesis  algorithms  and  show  the 
rougher  texturing  of  the  MPD  method  for  identical  fractal  dimensions.  The 
original  surfaces  were  constructed  using  floating-point  representation;  how¬ 
ever,  they  were  byte-scaled  for  3-D  plotting — ^thus  the  z-axis  range  of  0  to 
255  in  figures  6  to  9. 


Figure  1.  Plot  of  global 
mean  for  each  512  x 
512  block  processed 
with  PE  estimation 
compared  to  global 
fractal  dimension  used 
in  original  synthesis  of 
512  X  512  constructs. 
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Figure  2.  Synthetic 
fractal  surfaces  and 
corresponding  DFTs 
created  by  single 
realization  PSD  and 
MPD  methods  with 
fractal  dimension  2.0. 
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Figure  3.  Synthetic 
fractal  surfaces  and 
corresponding  DFTs 
created  by  single 
realization  PSD  and 
MPD  methods  with 
fractal  dimension  2.2. 
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Figure  4.  Synthetic 
fractal  surfaces  and 
corresponding  DFTs 
created  by  single 
realization  PSD  and 
MPD  methods  with 
fractal  dimension  2.4. 


Figure  5.  Synthetic 
fractal  surfaces  and 
corresponding  DFTs 
created  by  single 
realization  PSD  and 
MPD  methods  with 
fractal  dimension  2.6. 
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Figure  8.  Typical 
fractal  surfaces 
synthesized  by 
ensemble  average  PSD 
and  MPD  algorithms 
using  a  fractal 
dimension  of  2.4. 


Figure  9.  Typical 
fractal  surfaces 
synthesized  by 
ensemble  average  PSD 
and  MPD  algorithms 
using  a  fractal 
dimension  of  2.6. 


Table  4.  Global  mean, 
variance,  and  mamse 
for  each  512  x  512 
block  synthesized  by 
three  techniques  and 
processed  with  PE 
estimator.  (Expected 
fractal  dimension  is 
simply  value  used  in 
original  constructs.) 


Table  4  shows  the  average  values  for  the  estimates  in  tables  1  to  3,  and 
these  average  values  for  the  fractal  dimension  are  plotted  in  figure  1.  The 
key  in  figure  1  indicates  the  three  classes  of  images  produced  with  “Ex¬ 
pected  Fdim,”  indicating  a  perfect  estimate  based  on  the  fractal  dimension 
used  for  synthesis,  “PSD  ens.  Fdim,”  indicating  the  results  of  estimating  the 
PSD  ensemble  average  surfaces,  “PSD  sr.  Fdim,  ”  indicating  the  results  of 
estimating  the  PSD  single  realization  surfaces  and,  last,  “MD  point  Fdim,” 
containing  the  results  of  estimation  applied  to  the  MPD  synthesized  sur¬ 
faces.  These  results  show  quite  clearly  the  tracking  of  the  fractal  dimension 
by  the  simple  estimator  over  the  desired  range  with  an  increasing  deviation 
from  the  true  fractal  dimension  after  a  2.5  fractal  dimension.  Results  similar 
to  this  were  achieved  with  other  seed  values  used  in  the  random  number 
generator  for  producing  PSD  surfaces  as  well. 

The  algorithm  was  developed  using  the  PV-Wave  command  language  and 
is  currently  not  optimized  for  speed.  The  execution  time  for  the  algorithm 
has  not  been  determined.  The  unoptimized  PV-Wave  version  executes  an 
estimation  in  260  s,  including  the  system  overhead;  this  time  is  similar  to 
Moghaddam’s  [29]  optimized  C  version. 


Fractal  Expected 

surface  fractal 

type  dimension 

Calculated 

fractal 

dimension 

Variance 

Modified 
average  mean 
square  error 

MPD  method 

2.0 

2.1380 

0.0622 

0.0190 

2.2 

2.3333 

0.0558 

0.0178 

2.4 

2.6115 

0.0491 

0.0447 

2.6 

2.8884 

0.0428 

0.0831 

2.8 

3.1247 

0.0381 

0.1055 

3.0 

3.2968 

0.0344 

0.0881 

PSD  method  ensemble  average 
2.0 

1.9938 

0.0373 

0.0011 

2.2 

2.1059 

0.0330 

0.0111 

2.4 

2.3518 

0.0312 

0.0050 

2.6 

2.6542 

0.0301 

0.0054 

2.8 

2.9539 

0.0298 

0.0253 

3.0 

3.2145 

0.0309 

0.0463 

PSD  method  single  realization 

2.0 

2.0549 

0.0429 

0.0047 

2.2 

2.2119 

0.0403 

0.0025 

2.4 

2.4532 

0.0364 

0.0042 

2.6 

2.7281 

0.0341 

0.0172 

2.8 

2.9998 

0.0333 

0.0407 

3.0 

3.2391 

0.0335 

0.0579 

Values  are  the  averages  from  tables  1  to  3  with  the  calculated  fractal  dimensions 
plotted  on  plot  1. 
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4.  Conclusions 


The  PE  method  remains  a  possible  candidate  for  “real-time”  processing  of 
the  fractal  dimension  with  subsequent  texture  classification  and  segmenta¬ 
tion.  The  relatively  low  variance  and  mamse  demonstrate  the  algorithm’s 
accuracy  over  the  limited  range  of  interest. 

The  effect  of  the  rectangular  window  was  apparent  at  the  higher  dimensions 
due  to  side-lobe  leakage  and  therefore  must  be  questioned  for  the  segmen¬ 
tation  of  tree  clutter  [13]  and  probably  other  types  of  clutter  that  exhibit  a 
high  fractal  dimension.  To  circumvent  this  problem,  an  adaptive  estimation 
routine  may  be  useful  that  uses  the  results  of  Bovik  and  Super  [9],  together 
with  the  simple  algorithm  given  here. 

The  estimator  results  show  no  significant  difference  between  estimating  the 
ensemble  average  as  opposed  to  the  single  realizations.  This  should  be  ex¬ 
pected  since  the  DFTs  of  an  8  x  8  local  window  derived  from  the  ensemble 
average  are  indistinct  from  those  DFTs  taken  from  single  realizations.  The 
only  differences  I  observed  were  in  taking  global  estimates  (unpublished 
finding)  of  the  surfaces  and  the  slight  increase  in  the  variance.  This  indi¬ 
cates  the  viability  of  the  estimation  technique  to  “real  terrain”  (single  real¬ 
ization)  and  not  just  well-behaved  laboratory  models. 

The  presence  of  artifacts  within  a  mid-point  displacement  rendering  [5] 
probably  contributes  to  the  extra  power  distributed  along  the  spatial- 
frequency  axis  (see  fig.  1,  lower  right  comer),  thus  perturbing  the  estimation 
of  the  fractal  dimension.  A  path  to  further  resolving  this  issue  is  in  generat¬ 
ing  fractals  using  the  random  additions  refinement  to  the  MPD  method  and 
performing  the  estimation  [5]. 

Future  work  will  involve  the  optimization  of  the  PSD  synthesis  routine  and 
PE  estimation  algorithm  in  C  and  should  be  included  in  the  KHOROS 
image-processing  software  for  future  use.  An  optimized  version  of  the  PE 
algorithm  could  reduce  the  processing  time  by,  at  minimum,  a  factor  of 
four.  Accurate  measurements  of  the  optimized  code  will  be  performed,  as 
well  as  the  texture  analysis  and  segmentation  of  synthetic  composite  sur¬ 
faces,  natural  surfaces,  and  SAR  images. 
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